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Computer simulation of model cohesive powders: 
Plastic consolidation, structural changes and elasticity under isotropic loads. 
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The quasistatic behavior of a simple 2D model of a cohesive powder under isotropic loads is 
investigated by Discrete Element simulations. We ignore contact plasticity and focus on the effect 
of geometry and collective rearrangements on the material behavior. The loose packing states, as 
assembled and characterized in a previous numerical study [Gilabert, Roux and Castellanos, Phys. 
Rev. E 75, 011303 (2007)], are observed, under growing confining pressure P, to undergo important 
structural changes, while solid fraction $ irreversibly increases (typically, from 0.4-0.5 to 0.75- 
0.8). The system state goes through three stages, with different forms of the plastic consolidation 
curve, i.e., $ as a function of the growing reduced pressure P* = Pa/Fo, defined with adhesion 
force _fo and grain diameter a. In the low-confinement regime (I), the system undergoes negligible 
plastic compaction, and its structure is influenced by the assembling process. In regime II the 
material state is independent of initial conditions, and the void ratio varies linearly with log P [i. e. 
A(l/$) — AA(logP*)], as described in the engineering literature. Plasticity index A is reduced in 
the presence of a small rolling resistance (RR). In the last stage of compaction (III), "I> approaches 
an asymptotic, maximum solid fraction $max, as a power law, $max — ^ oc (P*)~", with a ~ 1, 
and properties of cohesionless granular packs are gradually retrieved. Under consolidation, while the 
range ^ of fractal density correlations decreases, force patterns reorganize from self-balanced clusters 
to force chains, with correlative evolutions of force distributions, and elastic moduli increase by a 
large amount. Plastic deformation events correspond to very small changes in the network topology, 
while the denser regions tend to move like rigid bodies. Elastic properties are dominated by the 
bending of thin junctions in loose systems. For growing RR those tend to reduce to particle chains, 
the folding of which, rather than tensile ruptures, controls plastic compaction. 

PACS numbers: 45.70.-n,81.40.Lm,61.43.Hv,83.10.Rs 



INTRODUCTION 



Cohesive granular materials are present in many nat- 
ural or industrial processes, the understanding of which 
requires studies of their rheology under small confining 
pressures, when tensile intergranular forces play a major 
role. In such cases cohesive materials exhibit specific fea- 
tures that do not exist in cohesionless grain assemblies, 
such as the ability to form stable structures at low den- 
sity and the sensitivity to stress intensity, as opposed to 
stress direction. Macroscopic constitutive laws and phe- 
nomenological tools have been developed and used in sev- 
eral engineering fields: mechanics of cohesive soils (clays 
and silts) [h,, i3, i^ iJ], metallic powder processing 3], 
modeling and treatment of ceramic powders [a, 0, B j 
handling of xerographic toners ll(|. One simple mate- 
rial is the assembly of wet beads [ll|, \1^ 13 1 , in which 
some microscopic observations are possible 12 UM ■ How- 
ever, wet grain packs are only slightly less dense than 
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dry ones, and do not enable the study of loose structures 
obtained with powders. In general, the behavior of ma- 
terials under proportional load (oedometric or isotropic 
compression) is characterized by the consolidation curve, 
which describes the irreversible compaction under grow- 
ing stress |l!]. Density can increase by factors of 3 or 4 
under growing load. 

Although numerical simulations have been widely used 
for several decades T4| to investigate microscopic mech- 
anisms and classify mechanical properties of granular 
systems, studies of cohesive materials are still far less 
common, and almost exclusively limited to dense ma- 
terials. Thus, the effects of capillary cohesion in wet 
sand or bead packs have been simulated [T^ [la], as 
well as the compaction of ceramic and metallic pow- 
ders [13, M, \M, M, M, M, [11 to states of very high 
density, or the behavior in shear tests of 2D dense cohe- 
sive packs with plastic deformation of contacts [2J, [2^ . 
Loose structures formed by particles packed under grav- 
ity and stabilized thanks to adhesion have been simu- 
lated [20] • Of particular relevance to the present study, 
among the very scarce numerical studies of loose pack- 
ings [23] stabilized by cohesion and of their collapsing 
under growing loads, arc the works by Bartels, Kadau, 
Wolf et al. [23, [23, [30, [sii on the oedometric compres- 
sion of granular assemblies with initial low densities. This 



research group studied a dynamical compression regime, 
and observed a shock wave propagating through the sam- 
ple. Shear flows of cohesive granular materials have also 
been simulated [H [H [H lH . 



II. MODEL MATERIAL AND SIMULATION 
PROCEDURES 

A. Definitions and basic equations 



In a previous article [36J, hereafter referred to as pa- 
per I, we studied by numerical simulation the assembling 
process, the structure and the force patterns of a model, 
two-dimensional (2D) cohesive granular material in loose 
equilibrium configurations. We now investigate the me- 
chanical behavior of the same model granular material 
in isotropic compression and pressure cycles, as well as 
the evolution of various characteristics of intermediate 
equilibrium states as plastic compaction proceeds. 

As in paper I, we keep the external pressure as the 
main control parameter. The adhesive strength Fq in 
contacts sets a force scale in the material behavior, and 
hence (in 2D) the reduced pressure, defined as 



P* 



aP 



(1) 



in which a is a typical grain diameter, is a crucial di- 
mensionless state parameter. The main objective of the 
present paper is the study of the process by which, as 
pressure is increased, cohesion-dominated loose struc- 
tures, for which P* ^ 1, get irreversibly compacted 
as P* increases until pressure dominates (P* ^ 1). 
Such a compaction was numerically observed e.g., in 
Ref. [31]. However, our approach produces homogeneous, 
isotropic, equilibrium configurations under varying load 
and is therefore apt to provide more detailed information 
about the connections between macroscopic constitutive 
laws and microstructural or micromechanical features. 

The present paper is self-contained and can be under- 
stood without reading paper I. A summarized description 
of the material properties and of the initial configura- 
tions (studied in paper I) is provided in Section [TTl The 
macroscopic material response in isotropic compression, 
with the possible influence of the initial state proper- 
ties, is studied in Section [IIII Then, various microscopic 
aspects of the consolidation process are investigated in 
the sequel: density correlations (with their fractal behav- 
ior over some length scale (361] ') are investigated in Sec- 
tion |lVl force networks and force distributions are dealt 
with in Section |Vl while Section IVII focuses on elastic 
moduli. Section IVIII discusses qualitatively some micro- 
scopic aspects of the consolidation behavior. The final 
section, part IVIIIi summarizes the results and suggests 
directions for future work. Sections IIVI and |V] can, at 
first, be read independently from each other. The same 
remark applies to Sections IVII and IVIII 



The material and the simulation method are identical 
to those of paper I [36| , which the reader might refer to for 
additional technical details, and for a physical discussion 
of some of the model ingredients. For the sake of com- 
pleteness, we however provide a summarized description 
below. The contact law is an elaboration of the often 
employed spring-dashpot model with Coulomb friction, 
in which two additional ingredients are introduced: an 
attractive force and, possibly, some resistance to rolling 
at contacts. The model material is a 2D assembly of 
disks with diameters uniformly distributed between a/2 
and a, enclosed in a rectangular cell with periodic bound- 
ary conditions in both directions. Both lengths Li, L2 
defining the cell size and shape are variable, and sat- 
isfy equations of motion designed to impose given values 
of diagonal stress components 0-1=0-2 = P- Stresses 
are controlled by a variant of the Parrinello-Rahman 
method [37|. In equilibrium, both diagonal stress com- 
ponents a on {oi = 1,2), with the convention that tensile 
stresses are negative, are given by the standard formula 
{A is the sample surface area): 



1 



l<i<j<N 



(2) 



In ([2]), the r.h.s. sum runs over all interacting pairs i,j 
among the N disks of the system, Fy is the force trans- 
mitted from grain i to its neighbor j and vector Tij points 
from the center of i to the center of j (with the suitable 
nearest image convention to account for periodicity) . The 
implementation of stress-controlled simulations is such 
that the cell length La along direction a increases or 
decreases if ctq is larger (respectively: smaller) than its 
prescribed value. 

As usual in molecular dynamics applied to granular 
materials (also known as the "discrete element method" ) 
particles have rigid body kinematics and their motion is 
governed by Newton's equations. 



B. Interaction law^ 

Grains interact with forces of elastic, adhesive, fric- 
tional and viscous origins. The static part of the normal 
component F^ of the force transmitted by grain i to its 
neighbor j is a function of hij, the distance separating 
disk perimeters. A negative hij means that the grains 
overlap, in which case they repel each other with a nor- 
mal elastic force F^'"* = —K^hij. This force vanishes 
whenever hij > 0. (Overlap hij < is, of course, a 
numerical representation of the physical contact deflec- 
tion). The repulsive elastic force is supplemented with an 
attractive term -F^'*"* , equal to —Fq for contacting disks 



{hij < 0). F^'^-' has a finite range Do, fixed to 10~^a, 
and varies linearly between —Fq and zero as hij grows 
from to Dq . Fq is the maximum tensile force a contact 
might support without breaking off. The normal contact 
law thus introduces a force scale, and a dimensionless pa- 
rameter, the stiffness parameter, k = gKn/Fq. k charac- 
terizes the amount of elastic deflection /ig under contact 
force Fq, relative to grain size a (ho/a = k~^). k is set 
to a large value, k = 10^, so that the elastic deflections 
in contacts remain so small that they can be neglected in 
comparison to all other length scales in the problem (in- 
cluding interstices between neighbors [3a])- The packing 
geometry can be regarded as that of an assembly of rigid 
grains (as formally dealt with in the "contact dynamics" 
simulation method used in [3l|). 

To the static contributions F^ and F^ to the normal 
force we add a viscous damping term opposing the rel- 
ative normal velocity of i and j when the disks touch 
(hij < 0), corresponding to a constant, positive normal 
coefficient of restitution cn in binary collisions if Fq is 
set to zero, cn is set to a low value, bn = 0.015 in 
our simulations. In the presence of attractive forces the 
apparent restitution coefficient in a collision will depend 
on the initial relative velocity. For small kinetic ener- 
gies the particles will eventually stick to each other. The 
minimum receding velocity for two particles of unit mass 
(the unit mass is chosen equal to the mass of a disk of 
diameter a) to separate is V*^/2, with 



the rolling friction threshold is not reached; and a rolling 
friction coefficient, fiR with the dimension of a length, 
setting the maximum absolute value of the rolling mo- 
ment Tji to (J-rF^ , proportional to the elastic part of the 
normal force. The implementation of this rolling law is 
analogous to that of the tangential one, with the rolling 
moment and the relative rotation respectively replacing 
the tangential force and the relative tangential displace- 
ment. A contact for which the total normal force is equal 
to zero in equilibrium, with F^ = K^ho — Fq, may 
transmit a rolling moment Tr with \Tr\ < (J-rF^. Since 
point contacts do not transmit torques, the rolling resis- 
tance stems from the irregularity of grain surface. Two 
contacting grains touch each other, in general, by two 
points (in 2D), which are separated by some microscopic 
distance / that is characteristic of the particle shape. fiR 
should be proportional to I, and Kr proportional to P. 
We set ^R — ^l and Kr = K^P, with, in most calcula- 
tions with RR, / = a/100. 

Table |T] summarizes the values of parameters used in 
most simulations, in dimensionless form. Some calcula- 



^^ 


ejv 


Kt 

'^ -77- 
Kn 


Do 
a 


Kr 

KNa^ 


a 


0.5 


0.015 


10^ 1 


10'^ 


10-^ 


or 0.005 



TABLE I: Values of dimensionless model parameters used in 
most simulations. 



V* = ^FqDq. 



(3) 



The elastic tangential force in contact i,j, F^ , is to be 
evaluated incrementally. In case of no tangential sliding, 
it varies linearly with the relative tangential displacement 
at the contact point, involving a tangential stiffness con- 
stant, Kt- In the case of sliding, which occurs when the 
elastic law would cause F^ to pass one of the Coulomb 



bounds ±^_F^'*-' 



then Fj! 



jY , uiicii ± J. stays equal to ±F^^-' . The 
relative tangential displacement at the contact point in- 
volves displacements of disk centers and rotations. The 
Coulomb condition introduces the friction coefficient, fi. 
It should be pointed out that it applies to the elastic 
repulsive part of the normal force only. Thus, a pair 
of contacting grains with hij equal to Fq/Kn = hQ, the 
equilibrium distance, such that the sum of elastic and ad- 
hesive terms vanishes, can transmit a tangential force Ft 
such that \Ft\ < /i-Fg- (The importance of this feature of 
the contact law for collective properties macroscopic be- 
havior of particle assemblies was stressed in paper I for 
isotropic, static states, and in Ref. [3J] in steady-state 
shear flows). All simulations reported here were carried 
out with fi = 0.5. 

We studied the influence of rolling resistance (RR) at 
contacts, which is modeled as in [39|. Two additional 
parameters are necessary: a rolling spring constant, Kr, 
with dimension of a moment, expressing proportional- 
ity between relative rotation and rolling moment (i. e., 
a torque concentrated at the contact point), as long as 



tions were also performed with larger RR (up to Z = a, 
Hr = 0.5a). 



C. Initial states 

In paper I, two extreme cases were studied in the as- 
sembling stage of cohesive packings under low P*. First, 
an A'^-particle sample of hard-disk fluid is prepared at 
solid fraction <&/ in a fixed cell. Then, in type 1 sys- 
tems, velocities are set to zero and the external pres- 
sure control is started, until an equilibrium is reached 
under F* = 0.01. The other procedure, by which type 2 
samples are prepared, is meant to represent the opposite 
situation, in which aggregation is much faster than com- 
pression. Thus, while the cell size is fixed and the solid 
fraction stays equal to $7, grains are attributed random 
(Maxwell-distributed) velocities and left to interact and 
aggregate until all N of them join to form one unique 
cluster. The system is then equilibrated at P* = 0, and 
compressed to P* = 0.01. To limit the influence of dy- 
namical effects, the strain rate is requested not to exceed 
a maximum value Cmax during compression. We express 
this condition with the natural inertial time associated 
with the characteristic force Fq: (m is the mass of a disk 
of diameter a) 



Tn 



(4) 



Sample type 


No cohesion 


Type 1 


Type 2 


N 


1400 


1400 


1400 5600 10976 


Number of samples 


4 


4 


5 3 1 


Lowest pressure 


P/Kn = 10"' 


P* = 0.01 


P* = 0.01 


$ (no RR) 


0.811 ±0.001 


0.723 ± 0.001 


0.472 ± 0.008 


$ (RR) 


0.805 ± 0.002 


0.688 ±0.001 


0.524 ± 0.008 



TABLE II: Set of granular samples used as initial equilibrated configurations in simulations of isotropic compression (with 
material parameters of Table HJ . 



defining a diniensionless inertia parameter 



Ia = e 



niaxJ^ 0- 



(5) 



la is set to 0.05 in our simulations. The main set of 
samples of types 1 and 2 (the latter coinciding with "se- 
ries A" in paper I), to which some non-cohesive ones 
are added for comparison, is listed in Table |lll in which 
the number of available configurations of different sizes 
is provided, along with solid fraction under the lowest 
nonzero pressure. All configurations are prepared both 
with {iir/o. = 0.005) and without {^a/a — 0) RR, with 
the parameters of Table HI The initial solid fraction is 
$/ = 0.36. Type 2 systems are also available under 
P* — 0, right at the end of the aggregation stage [3^ . 
but we regard this intermediate stage as part of the initial 
packing process and focus our study on higher pressures 
(as apparent in Table [Til the compression from zero pres- 
sure to P* — 0.01 involves a large density increase, and 
important changes of the microstructure are reported in 
paper I). Distant interactions between grain pairs sep- 
arated by a gap smaller than Dq are scarce, and "rat- 
tlers", i.e., isolated, free grains with no interactions, are 
absent in cohesive systems because of the initial aggre- 
gation process. Coordination numbers under P* = 0.01 
are typically z ~ 3.1 without RR, and z ~ 3.0 with RR, 
for both type 1 and type 2 cohesive samples. Additional 
details about those equilibrium configurations under low 
pressure can be found in paper I. 

The assembling stage of type 2 systems also depends on 
the initial velocities given to the grains before they form 
aggregates (the "granular temperature" of the original 
"granular gas"). The relevant dimensionless parameter 
is the ratio of the initial mean quadratic velocity Vq to 
the characteristic velocity V* defined in ([3]). Vq/V* is set 
to 9.5 for the main sample series of Table HIl The value of 
Vq/V* was shown in paper I to have a strong influence on 
the initial coordination number z at P* — in samples 
with RR: whereas z is larger than 3 for Vq/V* — 100, it 
approaches 2 for small Vb, of order y*/10, in which case 
the loopless structures of geometric ballistic aggregation 
models are retrieved. However, this effect is strongly re- 
duced after the compression step to P* = 0.01. 

In the following, unless otherwise specified, all results 
will pertain to the systems of Table HIl and measurements 
will be averaged over all available samples, error bars 
on graphs extending to one sample to sample standard 
deviation on each side of the mean value. 



D. Simulation procedures 

1. Equilibrium conditions 

One of the specificities of our simulations of cohesive 
packings under varying pressure is the approach, com- 
puting cost permitting, of the quasistatic material re- 
sponse, in which all configurations remain close to me- 
chanical equilibrium. Equilibrium conditions have to be 
stringent enough to enable an unambiguous identifica- 
tion of the force-carrying contact network and a study 
of its elastic properties. Due to the frequent occurrence 
of small contact force values, this requires forces to bal- 
ance with sufficient accuracy. We used similar criteria 
as in paper I, which, in a gree ment with other studies on 
cohesionless systems [33, |4Q|, were observed to provide 
adequately accurate force values. The tolerance levels on 
force and torque balance equations is expressed in terms 
of a typical intergranular force value Fi = max(i^Oj Pa)- 
A configuration is deemed equilibrated when (1) the net 
force on each disk is lower than 10~^Fi; (2) the total mo- 
ment on each disk is lower than 10~^Fia; (3) the differ- 
ence between imposed and measured stresses is less than 
10"^Fi/a; and (4) the kinetic energy per grain is less 
than 5xlO~*Fia. Those conditions being met, we could 
check that, in the absence of external perturbations (and 
of thermal motion) , no remaining slow motion, creep or 
aging phenomena were present in our systems: on waiting 
longer, only a very slow decrease of the remaining kinetic 
energy is observed. Furthermore, the computation of the 
stiffness (or "dynamical") matrix, see Sec. lIID"3l provides 
an additional stability check. 



2. Compression 

The sample series of Table [IT] are subjected to a step- 
wise compression cycle. In each compression step, ex- 
ternal reduced pressure P* is multiplied a constant fac- 
tor 10^/® ~ 1.334, and one waits until the new equilib- 
rium configuration is reached, with the criteria stated in 
Sec. HID II A condition of maximum strain rate is en- 
forced, in order to approach the quasistatic compression 
curve, as in the preparation process, on setting (see Eqs.IH] 
and 2]) la = 0.05. Parameter la, on replacing, in its def- 
inition, Fq by the force scale aP (in 2D) corresponding 



to the confining pressure is analogous to inertia param- 
eter / used to assess dynamical effects in steady shear 
flow [3J, [33] , or in the compression of non-cohesive gran- 
ular packings J38l . |4l| . The compression program is pur- 
sued until P* reaches the maximum value 13.33, above 
which negligible plastic collapse is observed. It should 
be noted that, thanks to the high value of stiffness pa- 
rameter K (see Sec. Ill A[) . the typical contact deflection 
qP/Kn at this highest pressure level is still very small. 
Then, the effect of decreasing P* back from its highest 
value to 0.01 is also simulated. As no large structural 
changes occur on decompressing the system, larger pres- 
sure jumps can be imposed on unloading. 

The simulations are computationally costly, as in some 
pressure steps equilibration times of order 100 To are re- 
quired, while the time step for the integration of the equa- 



tions of motion is a small fraction of yjm/ Km = TqJ \fK. 
This limits the size and the number of samples, and the 
use of small strain rates. Some tests of statistical signifl- 
cance and rate dependence of the results will be reported 
in Section [mi 



or vice-versa, one thus gets two separate measurements 
of the compliance matrix in our (statistically) isotropic 
systems, from which moduli G\\ and C\i are deduced, 
and hence the bulk modulus B = (Cn -f Ci2)/2 and the 
shear modulus G = (Cn — Ci2)/2. 



III. MATERIAL BEHAVIOR UNDER 
ISOTROPIC LOAD 



Compression and pressure cycle with 
non-cohesive material 



Non-cohesive systems of Table |TI1 initially obtained 
by isotropic compression of a granul ar g as (like the 3D 
sphere packings of e.g., Refs. [3a] and [44|), are subjected 
to a compression cycle, in which reduced pressure P/Kjq 
increases from its initial value Po/Kpj = 10~^, up to 
Pi/ Km — 1.33 X 10"^, and decreases back to 10^^. 

Typical results for the density of systems with and 
without RR are shown on Fig. [TJ Changes of solid frac- 



3. Computation of elastic moduli 

We observe that once samples are equilibrated accord- 
ing to the conditions of Section Hi D 11 then the Coulomb 
criterion \Ft\ < /ii^^r, as well as the rolling friction con- 
dition \Tji\ < (J^rF^ are satisfled as strict inequalities in 
all contacts. No contact is ready to yield in sliding, and 
with RR no contact is ready to yield in rolling either. 
This ensures that the response to small enough exter- 
nal load increments about a well-equilibrated state will 
be elastic and reversible. Elastic moduli express elastic 
response, i.e., with no effect of tangential or rotational 
sliding and no change in contact network topology and 
geometry. To compute elastic moduli, we build the stiff- 
ness matrix K of the contact structure (also taking into 
account the distant interactions) . K [36i | is a square ma- 
trix of order 3 A^ + 2 (the number of degrees of freedom in 
the system), depending on stiffness coefflcients K^ (re- 
placed by —Fq/Dq for the rare distant attractive bonds), 
Kt, Kr (with RR), and on network geometry. K is 
symmetric, positive deflnite (once the free translational 
motions of the whole sample as one rigid body are elim- 
inated) - and thus the stability of equilibrium states is 
checked. To compute elastic moduli, one solves a linear 
system of equations: 



K U 



(6) 



for the unknown displacement vector U, containing all 
particle displacements and rotations, as well as strains 
i^a)a=i 2- ^^^ right-hand-side of ([6]) contains external 
forces and torques applied to the grains, which are set to 
zero, and stress increments (A(Tq)^^j^ ^ (the same proce- 
dure is followed in ^2| with 2D disk packings and in [4^ 
with 3D sphere packings). On setting Acti = 1, A(T2 = 0, 



0.808 




0.806 



FIG. 1; (Color online) $ versus P/Kn in pressure cycle with 
1400 disk samples with and without RR. Blue dashed lines 
correspond to elastic response evaluated with the bulk mod- 
ulus from initial and highest pressure states. 

tion are very small (of order 10^"^, i.e., of order P/Km 
for the largest pressure), and nearly reversible (more 
than 90% of the density increase is recovered on de- 
compressing), as observed in Ref. [4l| with 3D sphere 
packings. The slight increase of bulk modulus as a func- 
tion of $ is due to the larger density of contacts under 
higher pressures. One typical feature of frictional, co- 
hesionless grain packs assembled by direct compression 
is the existence of a non-negligible population of "rat- 
tlers", i.e., particles that transmit no force (as observed 
e.g. in Ref. [H in 3D, or Ref. [12] in 2D systems). The 
fraction of rattlers xq thus exceeds 20% of the grains 
under Pg in systems with RR in the present case, and 
reaches 17% without RR. xq is reduced to 14% under 



P/Kn = 10~^. The backbone (force-carrying structure) 
is the set of non-rattler grains, characterized by coordi- 
nation number z* — z/(l — xq) [38|. z* increases with 
P, as rattlers get captured by the backbone and gaps 
separating neighboring grains close in compression. 

Changes of xq and z* are reversed on unloading (with 
some moderate hysteresis effect). The increase of z* as a 
function of P, above a minimum value Zg, which would 
correspond to P = 0, is sometimes described by a power 
law [45|. With such a fit we can estimate Zg, and we 
obtain values close to 3 with RR and about 3.12 without 
RR. z* varies by about 10% in the studied pressure in- 
terval. As in other simulations [SaEa, |43, the minimum 
coordination numbers stay above the "critical" value for 
rigidity, which is equal to 3 without RR and to 2 with 

RR M- 

Cohesionless systems under isotropic pressure cycles 
thus behave nearly elastically in an isotropic pressure cy- 
cle. As the pressure increases by more than 2 orders of 
magnitude, while remaining in the rigid limit of k ^ 1, 
only small and nearly reversible changes in density and 
in other internal state variables are observed, (see [41[ 
for a more detailed discussion). A small level of RR has 
little effect on density and material properties. 




FIG. 2: (Color online) Equilibrium configuration of a sample 
of 1400 disks with RR in initial state, under P* = 0.01, for 
which $ = 0.5132. Line thicknesses encode normal force in- 
tensities, red strokes depict compressive forces while tensile 
ones are colored in green, and forces equal to zero in blue. 



B. Compressing cohesive systems: general 
observations 



Once subjected to a pressure cycle, as specified in 
Sec HID 21 the material prepared in initially loose states 
(type 2 of Table |TT1) behaves as shown in Figs [2l [3] and g] 
As the pressure increases, so does the density, and the 
large pores present under low P* gradually disappear. 
The maximum packing fraction, ^„-,-^„ ~ 0.774 ± 0.001 
in that case, is quite reproducible. $max is smaller 
than the solid fraction of cohesionless systems (for which 
$ > 0.805, see Fig. HI). 

From the shape of <I>(P*) curves at growing P*, three 
regimes can be distinguished. At first, in a range of 
reduced pressure P* of the order of the first nonzero 
value (10~^), thereafter called regime I, $ remain ap- 
proximately constant: the contact network supports the 
growing pressure without rearranging. Then, in a sec- 
ond pressure interval which we shall refer to as regime II, 
a fast compression is observed. Density variations slow 
down in regime III, for P* of order unity, as a maximum 
solid fraction <f>„,ax is approached. On reducing the pres- 
sure, <& then remains very close to $„ax: the compaction 
is irreversible. 

The consolidation curve is similar to the ones obtained 
by numerical simulations in Refs. [2 3 . |3l|, on impos- 
ing uniaxial strains to loose packings prepared by an 
anisotropic ballistic aggregation process, although our 
study differs from these works in several respects (see 
Section |T|). Refs. @, [Ml focus on regime III, and on dy- 
namical compaction processes, with a shock wave prop- 
agating through the sample. The variations of solid 




FIG. 3: (Color online) Sample of Fig. [3 with $ = 0.6305, 
equilibrated under P* = 0.178 (different length and force 
units) . 



fraction $ versus P* are shown in Fig. [5l for three sam- 
ples of different sizes. Since all three curves are close to 
one another, we conclude that the macroscopic behav- 
ior is correctly captured in our simulations. Our results 
for $(P*) also resemble experimental curves obtained on 
different materials, such as metallic powders [5|, or xero- 




FIG. 4: (Color online) Same sample as on Figs. [5] and [31 
under the maximum pressure P* — 13.3. Solid fraction is 
$ = 0.7778. 



ing gradually form compact solids. For metal particles 
with d=10 /im diameter, one can estimate the pressure 
Fq/(P corresponding to P* = 1 to be in the 0.1 MPa 
range, so that the very large P* values in the compaction 
experiment reveal a different physical origin of density 
increase. The stiffness parameter, k, is also significantly 
smaller in such experiments, with the consequence that 
plastic phenomena cannot be ignored (for a definition and 
discussion of n in Hertzian sphere packings, see [Jdj). 
Contact plasticity dominates in the numerical studies of 
Martin et al. [Il [S [H [2l| , which focus on very high 
densities (beyond the random close packing value), when 
the material, due to sintering, turns into a porous com- 
pact. Hence only the early stages of metalpowder com- 
paction, in which densities are quite low [5| correspond 
to our simulations. In the case of the xerographic ton- 
ers studied in [T^, H, \^ , P* = 1 , as discussed in [3g| , 
rather correspond to P ~ 10 Pa. Nevertheless, the con- 
tact behavior, as investigated by atomic force microscopy, 
is likely to involve plastic effects [50, \Em, 112, ^M ■ 



C. Regime I: role of the initial assembling process 
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FIG. 5: Consolidation and decompression curves in 3 samples 
(with RR) with different numbers of grains, as indicated. 



graphic toner [10, Sa], at least in regimes I and II. Poquil- 
lon et al. [5|, in particular, in an experimental study of a 
metallic powder, explicitly distinguish three compaction 
regimes, with the material elastically resisting compres- 
sion in regime I, and then some plastic compaction, first 
attributed to particle rearrangement, as we observe, and 
later to contact plasticity. This latter effect, which is 
not included in our model, is likely to explain the dif- 
ference under high P* between many experiments and 
our results: experimental curves do not appear to ap- 
proach an asymptotic density, but witness ongoing com- 
paction up to the highest investigated pressure levels. 
In the case of metallic powders [5| , quite high pressures 
are applied (hundreds of MPa), and, as revealed by di- 
rect microscopic observations, particles fusing or sinter- 



As shown in [36| (paper I), and briefly recalled in 
Sec. IlICi assembling conditions have a considerable in- 
fluence on packing density and microstructure under low 
P* . It should be assessed to what extent those important 
differences in the initial configurations affect the plastic 
consolidation curve, and whether such a variability tends 
to disappear once the material undergoes significant com- 
paction. This issue is investigated in this section, in 
which the effects of various features of the preparation 
process are observed. The role of some micromechanical 
parameters is also discussed. 



1. Compaction and aggregation in ttie assembling stage 

The most important feature of the assembling process 
is the competition between compression and aggregation, 
which leads to the difference between systems of type 1 
and 2, as defined in [SO] and recalled in Section lll CI Type 
1 samples reach a considerably higher densities from the 
beginning, under low P*. Fig.[B]compares the subsequent 
consolidation curves. As type 1 systems are initially con- 
siderably denser, they are able to support larger pres- 
sures before rearranging, hence a wider regime I plateau. 
However, the pressure increase eventually reaches a high 
enough value to induce further compaction, and the con- 
solidation curve is then very close to that of type 2 sys- 
tems (the difference is actually smaller than the sam- 
ple to sample r.m.s. fluctuation). Within the accuracy 
and statistical uncertainty of our simulations, the differ- 
ence between initial states of types 1 and 2, although 
large, thus appears to disappear eventually upon plasti- 
cally compacting the material. 
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FIG. 6: Consolidation curve in type 1 and type 2 samples. 



FIG. 7: Consolidation curve with two different values of la 



2. Effects of first compression step and strain rate 

In paper I [3a] important changes between P* = 
and P* — 0.01 in type 2 configurations were reported, 
as solid fraction $ increases from <&/ ~ 0.36 to about 
0.5 (see Table HI]) . One way to limit the effects of this 
first compression step causes the most dramatic change 
is to reduce the strain rate, setting parameter la to a 
lower value. As shown on Fig. [71 displaying the con- 
solidation curve obtained in TV = 1400 systems with the 
usual value la = 0.05 and with the smaller one I^ — 0.01, 
lower inertial effects in the initial stage, while the equi- 
librium configuration at P* = 0.01 is prepared, result in 
a lower density and tends to turn the initial plateau of 
the <i>(P*) curve into a gentle ascending slope. Later 
on, as consolidation proceeds, very similar curves are 
obtained with both values of maximum dimensionless 
strain rate la (Fig. [7]), although the smaller error bars 
(representing sample to sample r.m.s. fiuctuations) wit- 
ness smoother changes and better reproducibility for the 
slower compression. It may thus be concluded that the 
quasistatic consolidation curve is quite reasonably ap- 
proached with the standard compression procedure de- 
tailed in Section IHni for which /„ = 0.05. 




FIG. 8: Consolidation curve: effect of initial agitation level 
in aggregation stage, and inffuence of RR parameter. 



low P*, with smaller coordination numbers. However, 
such a change in material properties does not only affect 
the initial, regime I part of the consolidation curve; it 
also alters the macroscopic mechanical behavior at larger 
densities: the slope of the consolidation curve is lower for 
larger RR. 



3. Effect of initial agitation and influence of RR 

The initial agitation velocity (or "granular temper- 
ature"), as expressed by ratio Vq/V* in the aggrega- 
tion stage strongly influences the coordination number. 
Figs. [5] and [5] show how this initial influence affects the 
beginning of consolidation curves and, once again, fades 
out later on. Consolidation curves are shown in Fig.[5]for 
two different values of Vq/V* , one tenfold as large as the 
standard value 9.5 used in the sample series of Table HIl 
and the other one smaller by a factor of 100. Fig.[9]shows 
the effect of Vq on coordination number. An increase of 
rolling resistance (with fin — 0.5 instead of 0.005), simi- 
larly to a decrease of Vb, stabilizes looser systems under 



4. Conclusion on initial states and regime I 

Fragile tenuous structures due to aggregation are easily 
perturbed and sensitive to many factors in low consoli- 
dation states. In general, all perturbations favor some 
kind of preconsolidation effect, inducing denser, better 
coordinated structures. These effects are reduced in each 
one of the following situations: (1) if one waits until large 
aggregates form before applying a conflning pressure; (2) 
if the initial agitation velocity Vb is decreased; (3) for 
slower compression processes, especially when the very 
first non-vanishing pressure value is imposed; (4) with 
larger RR levels. As the material is further compressed 
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FIG. 9: Same as Fig. [S] for coordination number z as a func- 
tion of P* . 



in (nearly) quasistatic conditions, the same macroscopic 
behavior is retrieved for given microscopic force laws [i.e., 
in cases (1) to (3)], irrespective of the initial perturba- 
tions affecting the beginning of the consolidation process. 
Though we did not vary the level of viscous dissipation 
in normal collisions, lower values are expected to induce 
larger inertial effects, similarly to a faster compression. 
On the other hand, viscous forces slowing down the mo- 
tion of grains relatively to a surrounding fluid (often an 
important physical effect in fine powders) could reduce 
the effects of the initial agitation. 

Regime I, with no plastic strain, is also observed in 
some experiments. For example, the response in uniaxial 
compression (i.e., ui > 0, (T2 = cs = 0) of loose aggre- 
gates of micrometer-sized silica beads assembled by bal- 
listic deposition - in that case, an anisotropic process in 
which particles are thrown onto a substrate - was studied 
by Blum and Schrapler [5J]. The deposit, with volume 
fraction $ ~ 0.15, resists a stress of 500 Pa before plastic 
compaction is observed, which corresponds to a "reduced 
stress", defined, in analogy with P*, as a\ = (Tiof/Fo of 



order 10~ . In the simulations of Wolf et al. [3i|J some 
finite initial pressure increment also has to be applied 
before plastic collapse is observed. 



D. Regimes II and III: 
intrinsic consolidation behavior 

Once the peculiarities of the sample preparation and 
first compression stage are erased, we refer to the mate- 
rial evolutionas the intrinsic consolidation behavior. In 
order to compare the shape of the consolidation curve to 
other observations more directly and quantitatively (and 
also for a more fundamental reason to be stated further) 
we subsequently describe it with 1/<I>, instead of $, as a 
function of log P* . This conforms to its traditional pre- 
sentation in the literature [ll, 13, 0, S Il0| > which often uses 
the void ratio, e = (l/*&) — 1. 



Once the regime I ends, we obtain linear variations of 
e or 1/$ with logP*: 



11 p* 



(7) 



where Pq and the corresponding solid fraction $0 are the 
coordinates of the point where the system behavior joins 
the intrinsic consolidation curve in the available samples. 
Parameter A, known as the plasticity index, is observed 
in our case to decrease as fiR increases from zero (Fig. [8|) . 
We have also observed that the value of this index is not 
affected by the friction coefficient: in that sense, /i just 
displaces the whole consolidation curve vertically [53|. 

As the maximum solid fraction ^^^^ is approached, 
Eq. ([7|) is no longer valid, and the asymptotic regime is 
better described with a power law, as in |3l| : 
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(8) 



with a constant A and an exponent a (close to 1 in our 
results). In order to describe the consolidation curve in 
regimes II and III with a unique functional form, we use 
the following relation: 



1 1 ^^ r"* 



1 — exp 



p* 



-| 1/a^ 



(9) 



which introduces additional parameters Pj* and a, and 
crosses over from Eq. ([T]), for P* ^ Pj*, to Eq. ([H), for 
P* » Pj*. Constant A in ([H]) is set to A/(2a) on using ^ 
for large P* values, and P* is directly related to ^^^^: 
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Fig. (TU] summarizes the definition and the role of all pa- 
rameters of relation ^. A fit of our data to relation ^ 



Region I: I/O ~ ftinction{prcparation) 
Region II: 1/*^, - X ln(P*/P*^) 

Region III: I/O + const (P*/P*)"" 
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FIG. 10: Schematic view of intrinsic consolidation curve with 
regimes II and III, and role of parameters introduced in 
Eq. (O. 
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FIG. 11: (Color online) Consolidation data and fit to Eq. Q, 
for systems with and without (small) RR. 



is shown in Fig. [TT] It should be noted that even a 
small level of rolling resistance changes the plasticity in- 
dex. Values of parameters are listed in Table IIIH where 
we also included the fit parameters for the sample with 
tJ-FtJO' = 0.5 corresponding to the data of Figs. [5] and [SI 



Mfl/a 


-Pq* '^O a <I>„ax Ct 





0.0237 0.469 0.349 ± 0.019 0.7808 0.91 ± 0.10 


0.005 


0.0316 0.515 0.194 ±0.004 0.7745 1.08 ±0.16 


0.5 


0.0178 0.382 0.25 ± 0.01 0.724 0.86 ± 0.24 



TABLE III: Values of parameters A, $„ax and a used to fit 
the consolidation curve in systems of Table Ull and in a sample 
with larger RR, with Eq. ((9}. Correspondingly, P* values are 
0.271 ± 0.033 without RR, 0.900 ± 0.064 for ^n/a = 0.005, 
and 2.6 ± 0.4 for ^lR,/a = 0.5. 

As the consolidation curve in region II, defined by pa- 
rameters A and Pq*, is observed not to depend on initial 
conditions, our simulations support the following inter- 
pretation: sooner or later in the process of quasistatic 
isotropic compression, the system joins, in the P* — <^ 
plane, a certain locus, corresponding to compressive plas- 
tic yielding. This locus, which acts as an attractor in 
isotropic compression, is a straight line on using coor- 
dinates InP* and 1/<I>. The value of Pq simply signals 
where, depending on the preparation process, the yield 
locus is reached. Table U gives the values of the parame- 
ters defining the intrinsic curve, and of pressure Pg* where 
it is first reached in type 2 systems of Table [III 

Consequently, in a system prepared at a lower density, 
it should be possible to observe a wider interval of the in- 
trinsic consolidation line. We could explicitly check this 
property in the case of one sample with N = 5600, for 
which the first nonzero equilibrium confining pressure in 
the loading history is equal to 2 x 10~^ instead of 10~^. 
This sample appears to have reached regime II sooner 



(around Pq — 10~^, or possibly below). The correspond- 
ing data points lie on the intrinsic consolidation curve 
(or, at least, within a distance smaller than error bars) 
identified on fitting the data of the main sample series, 
which had a larger first compression step (to P* = 10~^) 
and a larger value of Pq* (about 3 x 10^^). The yield 
locus can thus be extrapolated to lower pressures and 
densities, with the same plasticity index A. On assem- 




FIG. 12: Comparison of data obtained on the one low Pq 
sample (open triangles), and Eq. ^ (continuous line) with 
the parameters of Table UlTl as deduced from a fit of the data 
(black triangles) from the more systematic simulation series 
with larger Pg . 

bling cohesive aggregates with arbitrarily low densities, 
and on stabilizing them under very low initial pressures, 
it is conceivable (although increasingly difficult in nu- 
merical simulation because of the computational cost, as 
well as in experiments, because of the system sensitivity 
to perturbations) to create equilibrium structures with 
smaller and smaller densities and to explore an increas- 
ingly larger interval of the intrinsic consolidation curve 
in the limit of Pq — > 0. The corresponding solid fraction 
^Q would then also tend to zero. This limit is compatible 
with the functional form used in Eq. ([7]), while the use of 
the alternative form [10, E3I , 



p* 

^0 



4> — $0 = t/ In 



would lead to contradictions in the limit of Pq -^ 0. 



E. Unloading behavior 

On the $ versus P* curves we have been showing so 
far, that the unloading branch, down to P* — 0.01, shows 
very little density change. This property is actually sat- 
isfied on decreasing the pressure from other configura- 
tions in the compression process. Thus Fig. [I3| shows 



11 




# cycle 
1 


Path 
A 4B B' 


2 


A * C ^ C 


3 


A 4 D * D' 


4 


A ^ E * E' 


5 


C' ^ C # E' 





■ 








/ x^. 


0.7751 


- 








' ■ •- _ 


0.775 














■ 




y^ 


,i 


^'•'' 


0.7749 




V 


^^ 




- 



10 



p* 



FIG. 13: (Color online) Effect of different (isotropic) unload- 
ing/reloading histories on solid fraction. The direct consoli- 
dation curve with decompression from the highest pressure, 
as shown in previous sections, is ABCDEE' (path 4). On 
unloading along lines BB', CC, DD', the system does not 
rearrange. Such paths are reversible and do not alter the ma- 
terial state, since paths 4 (small black dots) and 5 (large, open 
pink circles) superimpose in P*, $ plane. 



FIG. 14: (Color online) Analog of Fig. [TJ for the unloading 
behavior of a sample with RR from P* = 13.3 to P* = 0.01. 
Dotted lines correspond to the elastic response of the highest 
pressure state and the final state (P* = 0.01). 



IV. CONSOLIDATION AND DENSITY 
CORRELATIONS 



that, if P* is reduced to the initial level 0.01 from differ- 
ent states on the consolidation curve, density changes are 
hardly noticeable, and $ stays very close to the maximum 
value reached at the largest imposed pressure P* in the 
past. Furthermore, it is checked (in the case of sequence 
4, drawn with open circles in Fig. [T3|) that the material 
might be reloaded, with no notable density change until 
pressure P* is reached. P* is known in soil mechanics as 
the consolidation pressure^ and a material in a state such 
that P* < P* is said to be overconsolidated. Upon in- 
creasing the pressure beyond the consolidation value P* , 
the density irreversibly increases, and this compaction 
is described by the same curve as in the absence of in- 
termediate pressure cycle: the recompression curve from 
C retraces back the same evolution from D to E. Thus 
the material behavior conforms to the plasticity of clays 
in isotropic compression |l|. All decompressing paths in 
the P*, $ plane, along which P* < P* , are reversible. 
More precisely, they are similar to the pressure cycles ap- 
plied to cohesionless systems (Fig. [T]), and they do not 
depart much from the linear elastic response, as shown 
on Fig. [TH 

For the largest P* values, adhesion forces are domi- 
nated by the confining stress and are nearly negligible: on 
setting Fq to zero in equilibrated systems under P* > 10, 
we could check that the granular assembly finds a new 
equilibrium configuration with very small displacements 
and hardly any change in the contact network. 



The gradual collapse of the initially open structure of 
loose systems, as visually apparent on Figs. [21 [Sj and [4] 
and witnessed by the consolidation curve studied in Sec- 
tion [iTTI can be characterized by the density correlation 
indicators introduced in paper I. 

The initial aggregation process was shown in paper 
I to result in a fractal structure of the density field 
over intermediate scales, between the grain diameter and 
some characteristic correlation length ^. In the presence 
of rolling resistance, even with the small value 0.005a 
adopted for /j,/j, the observed fractal dimension is com- 
patible with the result of the ballistic aggregation model, 
dp — 1.55. The ballistic aggregation model is purely 
geometric, and corresponds to the irreversible bonding 
of particles or aggregates in each collision, with contacts 
that are rigid in translation and rotation. This limit case, 
for which the coordination number is equal to 2, is ap- 
proached under low pressure [3a | with large RR or small 
Vo/V* . Better coordinated systems obtained with small 
RR and/or larger Vq/V* have the same fractal dimen- 
sion. Systems with no RR, on the other hand, are closer 
to dense objets with dp — 1-9 |36j . 

The limitation of the fractal behavior by an upper 
length scale ^ is a well-known geometric necessity in a 
large system with finite particle packing fraction <I>, be- 
cause (in 2D) a fractal structure of dimension dp < "2 
within a square cell of edge length L exhibits an appar- 
ent density proportional to L''^"^. In physically relevant 
circumstances, systems with a finite packing fraction <I> 
and a fractal structure over some distance range have a 
finite correlation length ^ above which the average value 
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of $ is observed. One then has <I> ex ^''^ ^ or 

^cx^'-i/^^-''^), (10) 

the prefactor being specific to the particular system stud- 
ied. Systems with size L 3> C can then be regarded as 
homogeneous packings of fractal "blobs" of (linear) size £,. 
Such ideas are quite generally used, and were applied to 
semi-dilute polymer solutions [531, to silica [5^] or poly- 
meric jSTJ gels, in computer simulations of aggregation 
models [5M|, and to various complex, supramolecular ob- 
jects like fat crystals [53| or asphaltene aggregates [60]. 

One may expect that the density increase caused by the 
collapse, under growing load, of the tenuous structures 
formed by cohesive packings corresponds to a decrease 
in the fractal blob size ^, while dimension dp still de- 
scribes the scaling of density correlation at smaller scale. 
One should then observe the scaling predicted in ([TOj). 
This implicitly assumes that the small scale structure of 
the packing is not affected by the compaction process, 
which essentially breaks long, thin junctions and fills the 
largest pores. A clue in favor of such a scenario is pro- 
vided by the results of Sec. IIII C[ which suggest that the 
same structure is obtained if the material is directly pre- 
pared with some value of $, or if it is assembled first in 
a looser state and then isotropically compressed, up to 
solid fraction $. 

To compute dp and ^, we measure the "scattering in- 
tensity" /(fc), i.e. the Fourier transform of the density 
autocorrelation function, as we briefly recall now (see pa- 
per I for more details). Density field x(r), taking values 
1 within particles and outside, is first discretized on a 
regular mesh, then Fourier transformed, thereby obtain- 
ing x(k). We then evaluate /(k) = |x(k)| /A, A being 
the cell surface area. Invoking isotropy, it is a function 
of fc = ||k|| alone. I{k) should then vary proportionally 
to fc~'^^ for a <C 27r/fc <C ^, and reach some plateau for 
k < 27r/C. 

This approach was used in paper I, and yielded the 
same fractal dimension, dp ~ 1.52 in systems with RR, 
under P* = (sohd fraction $/ = 0.36) and $ = 0.01 
(solid fraction $o — 0.524±0.008), while ^ decreased from 
£,1 = 9.3 ± 0.4 to Co = 5.1 ± 0.2. It should be noted that 
these values are roughly compatible with relation ([TOj) (as 
(6/^0)^"''^ = 1.4 ±0.1 is close to $o/$7 = 1.46 ±0.02). 

Fig. [T5| shows the scattering function for similar con- 
sohdation states shown in FigO {P* = 0.01), in Fig. [3| 
(P* — 0.178), and for P* — 1. These results are averaged 
over the four largest samples (with RR) of Table [TTl In 
spite of the error bars, I{k) exhibits the expected form, it 
is approximately constant below some crossover wavevec- 
tor 27r/C which increases with $, and then decreases, with 
slope —dp on a logarithmic plot. Pressure P* = 0.178 is 
the largest one for which this latter feature is clearly ob- 
served, and I{k) data corresponding to smaller pressures 
are intermediate between P* = 0.01 and P* = 0.178 
curves. The arrows on the plot signal the identified val- 
ues of wavevector 27r/C, which values have been estimated 
by means of the fit function for I{k) presented in paper I. 
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FIG. 15: Scattering intensity per unit area versus wave vector 
k. Results are averaged over the four largest samples (with 
RR) of Table ini 



The curve corresponding to P* = 1 - a flat, low scatter- 
ing signal - is typical of dense, homogeneous media with 
no fractal range for density correlations. 

In view of the small value of £ reached in the loos- 
est configurations (those with P* = studied in paper 
I), relation ITUj) is difficult to test from density correla- 
tion data. Another characteristic length scale for density 
inhomogeneities, used in paper I, is the (mass) averaged 
radius of gyration of pores. It may provide an alternative 
definition of a blob size £,', proportional to £,. We observed 
C ~ C at P* = 0.01, In fact, this equality works well 
under very low consolidations. However, under higher 
confining pressures we have observed that the definition 
of C gives lower values than £. Figure (THI is a plot of C 
as a function of pressure. 

Despite the restricted fractal range, our observations 
therefore confirm the validity of the "fractal blob" model, 
with a constant dp and a correlation length ^ decreas- 
ing as consolidation proceeds, until a final, homogeneous 
structure similar to that of cohesionless packings (albeit 
somewhat looser) is obtained. Other values of dp are 
likely to be observed with other assembling processes 
(such as e.g., diffusion- limited cluster aggregation). 

Values of C and dp do not, however, entirely determine 
the mechanical properties of the system. The response 
of an aggregate to some mechanical perturbation should 
depend on its connectivity which, as explicitly shown in 
paper I, is independent of its fractal dimension (systems 
with different iir and/or prepared with different values of 
Vq/V* have the same dp, but very different coordination 
numbers - see also Section [III C Sp . 

Results concerning blob sizes in systems without RR, 
for which dp ~ 1.9 |36|, are similar. Different stress 
states and mechanical conditions might also produce 
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FIG. 16: Average radius of gyration of pores, ^', versus P* . 



other types of loose structures. As an example, in 
the slow steady state shear flovi' of a very similar ma- 
terial simulated in [35| under normal reduced pressure 
P* = 1.25 X 10~^ and shear stress ai2 — 1.5P anisotropic 
structures with $ ~ 0.6 were observed. 



V. PROPERTIES OF EQUILIBRIUM FORCE 
NETWORKS 

A. Average normal force 

Formula ([2]), as explained in Ref. [15] and in paper I, 
leads to a simple relation between the average normal 
force (-FV) in equilibrium, pressure P, solid fraction $ 
and coordination number z: 



(F^ 



N 



TT{(f)P _ TiraP 
z$(d) ~ 9z$ ■ 



(11) 



We observed formula pT|) (involving the first and sec- 
ond moments of the diameter distribution) to be accurate 
in all simulated states despite some approximations in- 
volved [3g|. However, as stressed in paper I, relation (fTTj) 
fails to estimate the typical contact forces in the network 
under low P* . Those reach values of order Fq [iy| . Nor- 
mal forces of both signs (as visible on Fig [2] and Fig [3|) 
coexist and, to a large extent, compensate under low P*. 



B. Coordination numbers 

In initial low-pressure states, the coordination number, 
z, as shown on Fig. 1171 is nearly equally shared between 
the contribution z+ of compressive bonds and z_ of ten- 
sile bonds. A small population (zq per grain) of contacts 



carry forces equal to zero (within the numerical toler- 
ance for force equilibrium). Those contacts, in which the 
normal deflection h takes the equilibrium value ho for iso- 
lated pairs 36* , tend to be more numerous in the absence 
of applied stress, if the aggregation process avoids the 
building of hyperstatic (overbraced) structures. Their 
number is quickly reduced once aggregates made under 
P = are subjected to some external stress and start 
rearranging. 

The population of contacts loaded in compression in- 
creases along the consolidation curve until it dominates 
at large P*. Upon unloading, the initial proportion of 
tensile forces is first retrieved, and z_ is eventually, un- 
der low P*, larger than z+. 
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FIG. 17; Coordination numbers versus P* in compression 
cycle (a) without and (b) with RR. Both plots display, from 
top to bottom, z, z+, Z-, and zq. The error bars (not shown) 
are about the size of symbols. Arrows close to the curves 
indicate the compression branch of the pressure cycle. 

The total coordination number increases very little in 
the pressure cycle. Our observations thus contradict 
some statements in the literature [6l|, [62| relating z to 
$ (even in cohesionless systems, $ and z can vary inde- 
pendently [3g|). In the course of plastic collapse of loose 
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structures, as the solid fraction increases by more than 
50%, we observe the number of contacts to increase by 
5% in systems without RR, and by 12% with RR. Such a 
smah variation of z in plastic compression contrasts with 
the comparatively very fast change of z in the quasielas- 
tic compression of cohesionless packings, as observed in 
Section [ill Al in which z increases by more than 10% for 
minute density increases. 

As to the number of distant attractive interactions, i.e. 
pairs of neighboring grains separated by a gap smaller 
than Z?o (contributing to z_), it is initially very low (typ- 
ically 10 in a sample of 5600 particles), and then increases 
with P* but remains below 2% of the total number of in- 
teractions. 



C. Distribution of forces 

Normal force distributions are (roughly) symmetric 
about zero in initial states under low P* [ly|, as shown 
in Fig. [THl Under low P*, tangential forces of order Fq 
are also frequently observed [36|, and the angle between 
the total contact force F and the normal unit vector n is 
not constrained by the Coulomb condition, which applies 
to F -I- Fon rather than to F. This explains the typi- 
cal patterns of self-balanced contact forces in small grain 
clusters, where compressive and tensile forces of order Fq 
compensate locally, as might be observed on Fig. [51 The 
Coulomb condition applying to F, on the other hand, fa- 
vors alignments and "force chains". Self-stressed small 
clusters form spontaneously when the disks aggregate, 
except for large RR and/or smaU VnlV* ^^■ 

As consolidation proceeds, under growing P*, normal 
force distributions develop a wider positive (compressive) 
side (Fig. [T5)) . while the finite value for Fjv = — Po is 
characteristic of the failure of bonds in traction. Forces 
eventually scale proportionally to P* at large P* , like in 
cohesionless systems [4l|, |42l , as shown by Fig. [T9| When 
P* reaches values of several unities, the force distribution 
is similar to that of cohesionless packings, with an addi- 
tional dwindling population of tensile contacts (Fig. WJ\ . 
Force distributions in systems with small RR are quite 
similar to those shown in Figs. [T51 and [TOl 
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FIG. 18: (Color online) Probability distribution function 
P{Fn) of static normal force in contacts, versus Fm/Fq, in 
systems with no RR, for P* = 0.01 (black), P* = 0.178 (red), 
P* = 1 (blue), P* = 2.37 (green). Distribution widens as 
pressure increases as indicated by the arrow. P(Pjv) is also 
shown for P* = 0.01 for the overconsolidated state (OCS) at 
the end of the pressure cycle (pink dashed line) . 
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FIG. 19: (Color online) Positive wing of probability distribu- 
tion function of rescaled normal forces, Fn/P*, in systems 
with no RR, under P* = 2.37 (black crosses), P* — 5.62 (red 
square dots), P* = 13.3 (blue triangles). 



D. Forces in dense, overconsolidated states 

Upon decompressing to low pressure levels, some larger 
compressive forces {Fn/Fq reaching 2 or 3) survive and 
the distribution is not symmetric (Fig. \TE\i . Such effects 
of overconsolidation on contact forces are considerably 
larger than in cohesionless granular materials [4l| . As in 
the case of cohesionless systems |4l|, we observed that 
the decompression process tends to be affected by dy- 
namical effects if it is too fast, and the overconsolida- 
tion effects on force distributions tend to be erased if too 
many contacts open in transient stages. The results per- 
taining to overconsolidated states shown in Figs. [TSl [T71 



and [20] were obtained on simply reversing the stepwise 
compression program with the parameters indicated in 
Section [lID 21 (i.e. with as many steps in decompression 
as in compression). 

This final force distribution is similar to the one re- 
ported by Richefeu et al. 16] in simulations of packings of 
wet spherical beads, in which cohesion is due to capillary 
forces. After assembling the packing under a finite pres- 
sure and then decompressing to P = 0, these authors ob- 
serve that the particles tend to form small domains with 
only compressive or only tensile forces. Fig. [20] reveals 
quite similar patterns in overconsolidated states under 
P* = 0.01, with some predominance of the regions under 
tension, while compressive forces tend to organize more 
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FIG. 20: (Color online) Dense overconsolidated state of a 
sample under P* = 0.01 (with RR) at the end of the pressure 
cycle. Color code as on Fig. [2] with distant attractive forces 
(for which < h < Do) in blue. 



tion 2^22 of the contacts that join 2-coordinated disks. 
Such contacts are impossible in an equihbrium structure 
without RR. a;22 reaches 12% in large RR systems under 
low pressure (for $ in the 0.4 to 0.5 range), down to 1- 
2% in the main sample series of Table HIl with small RR 
{fj.R/a = 0.005). Thin, rigid strands of 2-coordinated 
disks might, however, be decorated by a side arm act- 
ing as a dead end for force transmission, and their me- 
chanical role is thus only partially captured on simply 
recording fraction X22. In the limit of z — > 2, which is 
approached under low pressure for large RR and/or low 
velocity Vq in the assembling stage, the force network has 
a vanishing number of loops and approaches isostaticity, 
as discussed in paper I [36|. Consequently, as compared 
to the case of small or no RR, systems with large RR 
under P* <C 1 exhibit narrower force distributions. For 
P* of order 10^^, normal forces above -Fb/5 or below 
— i^o/10 &re extremely scarce (with probability distribu- 
tion function P{F]\f) in the 10~^ range). Furthermore, 
with P* ~ 1, while compressive normal forces of order 
Fq are frequently observed, P{F]\[) remains below 10~^ 
for Fn -^ —Fq. This contrasts with the results shown on 
Fig. [13 the proportion of contacts on the verge of ten- 
sile rupture is much smaller in systems with large rolling 
resistance. 



often in strong force chains. Tensile contacts are more 
numerous than compressive ones after the pressure cycle 

(Fig.ini). 

To what extent overconsolidation effects on inner 
states influence the mechanical properties of cohesive 
granular materials (e.g., their response to shear stress) 
would deserve to be investigated. 



E. Effect of a large rolling resistance 

As reported in the previous paragraphs, the small level 
of rolling resistance used in most simulations reported 
here [nn/a = 0.005) has no large effect on force distribu- 
tions or force patterns. Yet, such a small RR significantly 
affects plasticity index A (see Table HIH) and changes frac- 
tal dimension dp (Section lIVp . 

In order to understand the mechanisms by which RR 
affects macroscopic behavior and geometry, we investi- 
gated the effects of a large RR {^n/a — 0.5) in a few 
1400-disk configurations. Rolling resistance favors force 
transmission along thin strands of particles, each of them 
in contact with two neighbors (such structures are shown 
in paper 1 [361, Fig. 20]). Single particle chains are eas- 
ily disrupted if ^in/a is small, but are quite frequent for 
such RR levels. While single particle chains are easily dis- 
rupted if pLRJa is small, they become much more frequent 
for large rolling resistance Thus coordination numbers 
may approach 2 (see Fig. [H]). The density and the length 
of such particle chains are also witnessed by the propor- 



VI. ELASTIC MODULI 

Elastic moduli are used in experiments [63| and com- 
puter simulations [43, [M] to express the response of gran- 
ular materials to small load increments. Their measure- 
ment, or that of wave velocities, is a non-destructive 
probe of the packing structure. Thus, in the case of co- 
hesionless bead packings, the simulations of [ijl showed 
that the moduli are sensitive to coordination number, 
which can vary independently of the solid fraction, and 
escapes direct observations [3g|. In the present case of 
possibly loose and poorly connected cohesive systems, 
those moduli also approximately describe the parts of 
the compression curves with no packing rearrangement 
(Fig. [H)l . like in cohesionless systems (Fig. [1]). 



A. Elastic moduli of cohesionless packings 

We first quickly describe the variations of elastic mod- 
uli in the cohesionless systems of Table |IT1 and their rela- 
tions to microstructural or micromechanical parameters. 
Fig-inHis a plot of bulk and shear moduli versus pressure. 
Values of moduli are very similar in systems without and 
with RR, and vary very slowly with ^r in the latter case. 
Unlike with Hertzian contacts, local stiffness constants 
-ftTjv, Kt do not depend on forces. Consequently, the in- 
crease of moduli with pressure is moderate. The results 
of Fig. [nH are typical of cohesionle ss g ranular systems 
with small coordination number |42l . |43| . |47[ . The evolu- 
tion of bulk modulus is correctly described by the simple 
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FIG. 21: (Color online) Bulk and shear moduli of cohesion- 
less systems (with RR), versus pressure in compression cy- 
cle. Voigt and Reuss bounds are shown as (red) triangles and 
(blue) round dots, respectively. Asterisks show values of B 
obtained on taking a larger rolling stiffness, Kr — 10~^KNa^ 
instead of Kn — 10~''' Kncl^, with the same contact network. 



estimation formulae recalled below in Sec. IVIBl and it is 
explained by the increase of coordination number. Shear 
modulus G, on the other hand, is somewhat anomalously 
low, witnessing the propensity of a rather poorly con- 
nected contact network {z* ~ 3.1 under P/Km = 10~^, 
without RR) to rearrange under small stress increments, 
if those are not proportional to the preexisting stresses. 

The evolution of elastic moduli in the unloading part of 
the pressure cycle (not shown on the figures, for clarity) 
very nearly reverses the effect of the first compression. 



B. Simple estimation formulae 

Bulk and shear moduli are traditionally estimated by 
the Voigt or mean field formula [4Jj ISa] , which gives up- 
per bounds [4^ B^ , G^ in terms of contact stiffness 
constants and coordination number z, based on the as- 
sumption that particle centers move like points of a ho- 
mogeneously strained continuum. In the present case one 
has: 



B' 



G^ = 



z$ [{(f) -f (d)2] Kn _ ^bz<^K_ 



N 



K 



N 



2Kn 



Att{(P) 



1127r 



(12) 



On deriving (|12p . similar approximations are used as 
for pTjl . The formulae are identical for systems with 
or without RR, and since we chose Kt — Kn one has 
also G^ = B^. 

For the bulk modulus, one may also write down a lower 
bound B^, the Reuss estimate l43|, based on the eval- 
uation of the elastic energy with trial forces in a load 



increment. The formula for B^ involves moments of the 
contact force distribution, specifically the following ratio: 



{F, 



Kt T 



Kt, 
Kn 



,r2) 



{F. 



N) 



(13) 



in which averages are taken over all contacts carrying 
static normal force Fn, tangential force Ft and rolling 
moment F (to be set to zero in the absence of RR). Us- 
ing (|lip . one has 



B' 



z^id^KN _ 27z^Kn 
2ii{(P)Z2 ~ 567rZ2 



(14) 



This approximation of the bulk modulus becomes exact 
when the force increments caused by an isotropic pressure 
increase are proportional to the preexisting forces [43| . 
and hence it tends to be accurate in systems with small 
degrees of force indeterminacy. The ratio of upper to 
lower bounds for B given by Eqs. p^ and (HU is 
55Z2/54, and the bulk modulus is therefore especially 
well predicted when the force distribution is not too 
wide [431, and ratio Z2 stays close to 1. Thus bulk moduli 
are rather successfully estimated (see Fig. I^Tj) by B^ or 
B^ in the cohesionless case of Sections IIII Al and IVI Al 
Force distributions have often been studied in cohesion- 
less systems, in which they are strongly constrained by 
the no-tension condition, and Z2 cannot reach large val- 
ues {Z2 < 1.5 in the present case). 



C. Elastic moduli in cohesive packings 

Elastic moduli as functions of P* during consolidation 
of cohesive systems are plotted in Fig. [511 Note the log- 
arithmic scale used for elastic moduli (unlike in Fig. BTjl . 
Both bulk and shear moduli are very low at small P*, 
which cannot be simply explained by the factor z$ ap- 
pearing in estimates P^ and P^ (z values, see Fig. [T71 
are similar to those of cohesionless systems while $ is 
twice as small at most). Those anomalously low moduli 
witness the propensity of the system to rearrange under 
isotropic as well as under deviatoric stress increments. 
Moduli in samples with RR (Fig. [HJs) have very similar 
values as in the absence of RR, although this may be 
partly coincidental, since they are quite sensitive to the 
value of rolling stiffness K^. 

On decompressing, the moduli (not shown in Fig. [22)) 
stay close to the value reached at the highest pressure. 

Mean field estimates B^ and G^ are both too large by 
factors of 30 to 50 in loose states. From ([TT|) the average 
normal force (Fn) vanishes as P* tends to zero, while the 
second moment is of order Fg . Moreover, as tangential 
forces are not limited by condition \Ft\ < /J^Fn, but by 
\Ft\ < IJ'iFN + Fo) instead, their contribution to the elas- 
tic energy is important (and so is that of rolling moments 
in systems with RR) . Coefficients Z2 thus reach values of 
order 10^ or 10'^ under low pressure, whence B^ / B^ ^ 1 
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FIG. 22; (Color online) Bulk and shear moduli of cohesive 
systems (a) without and (b) with RR, versus (growing) pres- 
sure. Same symbols and colors as in Fig. 1211 



which is impossible in cohcsionless systems. The Reuss 
bound for B is first (in regime I) too small by a large 
factor. Then, it seems to capture the evolution of the 
bulk modulus in regimes II and III of the consolidation 
behavior. Ratio B/B^ is reduced to about 2 for P* of 
order 0.1, and slightly decreases as compression proceeds. 
It should be recalled, though, that the Reuss formula es- 
sentially relates the bulk modulus to another unknown 
quantity, Z2. 



D. Elastic moduli and force indeterminacy 

The low value of the shear modulus in poorly coordi- 
nated cohcsionless packings under isotropic stresses (see 
Fig. m]) has been observed [43, S^l and argued [6a| to 
stem from its tendency to vary proportionally to the de- 
gree of force indeterminacy per unit area (or volume in 
3D) when it is small. As the latter (without RR) is pro- 
portional to (z* — 3)$(1 — xq), one should have 
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FIG. 23: (Color online) Elastic moduli (no RR) divided by 
A']v$(l ~ xq), versus z* , the coordination number without 
rattlers. Data with error bars, fitted with the dashed straight 
line, correspond to G in cohesionless systems. G and B in the 
cohesive material are respectively shown as (red) crosses and 
asterisks. 



Fig. [53] shows our cohesionless packings to abide by this 
law, as the linear variation of G* with z* would predict, 
within uncertainties, its vanishing for z* = 3. However, 
it is also obvious from Fig. [52] that the anomalous be- 
havior of both moduli in loose, cohesive grain assemblies 
are not simply explained by their low coordination num- 
ber, except perhaps for the shear moduli of the densest 
configurations (rightmost data points), which, after suffi- 
cient plastic compaction, become similar to cohesionless 
packings. A coordination number z* just above 3 (with- 
out RR) characterizes a "barely rigid" contact network, 
but such a global, average quantity does not account for 
the specific heterogeneities of loose cohesive packings. 



E. Contact forces in a small pressure increment 

At the microscopic level the elastic response to a small 
pressure increment AP determines contact force incre- 
ments as visualized in Fig. 1241 Very strong compres- 
sive force chains appear, while large parts of the system 
carry very small forces. On sorting the contacts by de- 
creasing contribution to the elastic energy of the force 
increments balancing AP, less than half of them (46%) 
contribute 95% of the energy. This proportion increases 
to about 65% in the densest configurations, to be com- 
pared to 68-70% in cohesionless systems. The configu- 
ration of Fig. [M] in a system with RR, has quite a few 
dead ends, i.e., sets of grains that are connected to the 
rest of the structure but do not belong to any percolating 
loop for force (or current) transport through the whole 
periodic cell. With RR, the force-carrying structure coin- 
cides with the backbone in the sense of ordinary (scalar) 
percolation theory. The force patterns of Fig. [M] dif- 
fer from those of Fig. [2] in which the equilibrium forces, 
prior to the application of AP are shown: some regions, 
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FIG. 24: (Color online) Force increments associated with elas- 
tic response in isotropic compression of system of Fig. [2] Con- 
tacts are ordered by decreasing contribution to elastic energy, 
and only the first 46% contact forces corresponding to 95% of 
the energy are drawn (colors as on Figs.[2]to|4ll. 



fuzzy ^-sized objects, the blobs. To discuss elastic prop- 
erties, the system is better represented as a network of 
"superbonds" of length ^, or effective beams (with which 
the elongated structures carrying stress in Fig. [Ml could 
be identified). In such a network, the dominant defor- 
mation mode is beam bending. The transverse deflection 
6 in bending of a beam of length ^, caused by a force 
F, is proportional to £^^F. Macroscopically, strains are 
of order e — S/£^, while F corresponds to stress a by 
F (X a£, in 2D. Consequently, the scaling of elastic mod- 
uli a/e with length ^ should involve a factor ^~'^. (Some 
possible corrections to exponent 3 are possible, although 
the appropriate value in, e.g., the case of percolation net- 
works of beams is very close to 3 [67|, 163] ) . As ^ varies by 
a factor of 3 or 4 within the scaling range (see Fig. [TB)l , 
relation B oc ^~^ would predict an increase of moduli by 
a factor of a few tens. 

Although this can be regarded as a fair estimate (see 
Fig. E^ . it should be admitted that the fractal range is 
very likely too restricted for such scaling laws to apply 
without important corrections. With sufficiently large 
rolling resistance, the "beams" can be reduced to single 
particle chains, which, as we now show, enables simpler 
analyses of their bending stiffness. 



G. The case of a large rolling resistance 



especially the isolated, self-stressed clusters where com- 
pressions and tensions of order Fq equilibrate, carry large 
forces but are bypassed in the transmission of the pres- 
sure increment AP. Dead ends contain "islands" of self- 
balanced forces resulting from the aggregation process, 
as directly visible on Fig. O but they do not participate 
in the transmission of stress increments and they do not 
contribute to elastic moduli. 

As consolidation proceeds, the repeated application of 
pressure increments clearly favors force chains over local- 
ized self-stressed clusters, and the force pattern adapts to 
the external pressure. Hence a closer similarity between 
the spatial distribution of equilibrium forces under pres- 
sure P and that of force increments caused by a small 
compression step AP, and a better performance of the 
Reuss estimate. 



F. Scaling with fractal blob size 

The inability of the approaches used in cohesionless 
systems to predict the elastic moduli of loose cohesive 
packings can be attributed to their ignoring the peculiar 
network geometry, which is the origin of the strong force 
concentration shown in Fig. [24l 

In view of the results of Section IIVI it is tempting to 
relate the elastic moduli to the variations of blob size ^. 
In scaling arguments about the density, the system can 
be regarded as a densely packed assembly of somewhat 



With large RR the prevalence of particle strands as 
force-transmitting structures (Sec. IV Ep influences elas- 
tic properties. As noted above, linear structures tend to 
deform like bending beams, with a compliance propor- 
tional to the third power of their length. In the case of 
single linear strands, connections with contact properties 
are easily made more explicit. Consider, e.g., a straight, 
linear chain of n identical disks of radius R, with n — 1 
contacts characterized by stiffness constants K^, Kt and 
Kn. Then, in the elastic regime, all intermediate disks 
can be suppressed and the interaction between the ex- 
treme ones, numbers 1 and n along the chain, can be 
replaced by an effective one between two disks of radius 



{n — l)R, and compliances 1/Kj^ , 1/Kj^ , 



and l/K^R^ 



for normal, tangential and rolling relative motion, with: 
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For large n the tangential compliance is much larger than 
the longitudinal and rolling ones, so that long chains be- 
have as beams, which essentially deform in bending. The 
local bending stiffness EI of the beam (i.e. the product 
of the material Young modulus by the moment of inertia 
of the beam section) corresponding to the chain of parti- 
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cles in the continuous limit is EI = 2RKfj. (This coeffi- 
cient expresses the proportionahty of bending moment to 
rotation angle gradient). For n ^ 1, the bending spring 
constant 3EI/P (expressing the transverse force to trans- 
verse deflection relationship) is correctly identified from 

Kj^' given in p6)) . using the length / = 2(n — 1)R of the 
straight n-particle strand. 

Remarkably, the bending elasticity of small linear 
strands of micrometer-sized colloidal particles bound by 
adhesive forces has recently been measured by means of 
optical tweezers [69]. Colloidal gels of polymer parti- 
cles [70, [Zii 12] should thus be modeled as cohesive par- 
ticle assemblies with rather large RR level. 

It is easy to check (consider e.g., two such chains join- 
ing at their ends at some angle) that for all strand shapes 
other than straight lines, the extremities will be coupled 
by spring constants of order K!^ for both longitudinal 
(parallel to end-to-end vector) and transverse relative dis- 
placements. Consequently, the macroscopic elastic mod- 
uli should be proportional to rolling stiffness constant 
Kif. Fig. [21] shows that this proportionality is approx- 
imately satisfied in the loosest states of a system with 
rolling friction fiR/a = 0.05, in which three different val- 
ues of Kji were used to evaluate the elastic response. 
Elastic moduli of denser states, however, depart from 




FIG. 25: (Color online) Bulk (filled symbols) and shear (open 
symbols) moduli, normalized by Kr, in low pressure states 
of a sample for which /iR/a = 0.05 and Kr = 10~ Knu 
(black squares). Results obtained on evaluating moduli with 
Kb, = W'^Kno^ and with Kn — W'^'^Kncl^ are respectively 
shown as red triangles and blue circles. 



this behavior. Therefore, the scaling of elastic moduli 
with typical strand length (as suggested in Section [VlF[) 
is limited to low consolidation states. With small or van- 
ishing RR, single particle strands are replaced by thicker 
junctions, which further restricts the consolidation pres- 
sure range for which elasticity is dominated by beam 
bending. 



VII. PLASTIC CONSOLIDATION 
MECHANISM: QUALITATIVE ASPECTS 

Cohesionless granular assemblies, if subjected to stress 
increments that are not proportional to initial stresses, 
essentially deform because the contact network gets re- 
peatedly broken and repaired [io, [z3- Macroscopic 
strains, once they exceed the very small scales associated 
with the response of given contact networks [40, |43| , thus 
result from a sequence of rearrangement events or micro- 
scopic instabilities, during which the granular packing 
loses its coherence and gains some finite amount of kinetic 
energy, even for arbitrarily slow applied stress changes. 
Collisions and appearance of new contacts stabilize the 
packing at the end of each microscopic rearranging event. 
This process gradually changes the topology of the con- 
tact network, and produces specific evolution of its fabric 
(orientation anisotropy). 

The mechanism of plastic collapse in isotropic compres- 
sion of loose cohesive assemblies with small or vanishing 
RR in contacts, as observed in the present study, is sim- 
ilar. Just like in cohesionless systems under shear [40], 
we expect the frequency of occurrence of rearrangements, 
along the loading path, to increase, and the correspond- 
ing strain jumps to decrease, as the size of samples grows, 
and thus the consolidation curve should be smooth in 
the thermodynamic limit. Due to the specific geome- 
try of loose systems, in which dense zones are weakly 
connected through thin arms, better connected, solid- 
like regions tend to move like rigid bodies, while frag- 
ile junctions break and rearrange, so that initially large 
holes gradually fill up. Fig. [26] illustrates this scenario. 
Displacements are depicted as arrows, pointing from the 
current positions to the ones reached in the next equilib- 
rium configuration in the stepwise compression sequence. 
The more densely packed, nearly rigid regions (marked 
with dotted lines) are easily identified by direct visual 
inspection. Fig. [^D also shows that the contact network 
undergoes relatively small topological changes, as more 
than 90% of contacts are conserved. The rate of con- 
tact change, and the evolution of coordination number 
with strain are significantly smaller than in cohesionless 
systems undergoing, e.g., shear deformation. During the 
compaction of loose samples the dense regions collide and 
slide past one another, along thin sheared zones where 
most of the broken contacts are found. 

In the case of large RR, the peculiar microstructure 
involving single particle chains might lead to a differ- 
ent deformation mechanism. Unlike multiply connected 
junctions, simple strands can yield in bending without 
breaking: they fold at some contact, where the rolling 
friction threshold is reached, thereby releasing bending 
elasticity. This mechanism is observed in experiments 
on single chains of colloidal particles [60, IZll ■ One thus 
expects fewer contact losses in plastic compression. 

To follow more closely the rearrangement sequences 
in the course of compaction, it is appropriate to moni- 
tor changes in the list of contacts during the motion be- 
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FIG. 26: Equilibrium particle positions in 1400 disk sample 
with small RR under P* = 0.032. Particle displacements to 
new configuration equilibrated under P* = 0.042 are shown 
as arrows (global density change A<1> — 0.05). Neighbor pairs 
for which contact opens are filled in grey. All other contacts 
(thin solid lines) are maintained. Dense regions moving ap- 
proximately like rigid solids are circled within dotted lines. 
Most lost contacts are situated near the boundaries of such 
solid-like particle lumps. 



tween two equilibrium configurations. As an example, let 
us consider the evolution between equilibrated states as 
P* increases from 0.177 to 0.237, and compare two sam- 
ples, one with small {nR/a = 0.005) and the other with 
large {^R/a = 0.5) RR. Table IIVI gives the changes in 
solid fraction and coordination number, and numbers of 
maintained, destroyed and created contacts in this com- 
pression step. Successive configurations separated by a 
fixed time interval At — O.lSTg are compared and Fig. [771 
plots the number of destroyed and created contacts as 
functions of time. For the same strain increment, con- 
tact losses, as a function of global strain, are significantly 
less frequent in the sample with large RR. This fact is re- 
fiected both in the data of Table HVl where global changes 
are recorded, between the initial and final states, and in 
those of Fig. [27l where successive changes over time inter- 
vals Ai are detailed. As a consequence, while the coor- 
dination number hardly changes during consolidation in 
systems with small or vanishing RR (see Fig. \Vf\ , it grad- 
ually increases from an initial value close to 2 to nearly 3 
in systems with large RR (Fig.IH]). The lesser importance 
of tensile contact rupture in the plastic compression of as- 
semblies with large RR is also witnessed by the normal 
force distribution (Section |VE|: forces approaching — Fq 
are quite scarce, as opposed to the situation in samples 



0.005 3.2 0.14 2084 (94.9 %) 112 (5.1 %) 115 (5.2 %) 
0.5 3.1 1.2 1679 (98.5 %) 26 (1.5 %) 46 (2.7 %) 

TABLE IV: Relative changes of solid fraction, A$, and of co- 
ordination number (Az), and numbers of maintained [N^^'), 
destroyed {N^~') and created (A'''^') contacts in a 1400 disks 
sample, with small or large RR, in the compression step be- 
tween P* = 0.177 and P* = 0.237. 
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FIG. 27: (Color online) Evolution of the contact number 
as a function of relative density increase. In sample with 
liRJo. = 0.005 the proportions a;+ and X- of gained and of lost 
contacts with respect to the previous recorded list are respec- 
tively shown with red square dots and triangles - the latter 
being connected with a dashed line. A similar code is used for 
a;+ and X- values in a sample with large RR {^R/a — 0.5), 
but with open dots, and in black. 



without RR (Fig. [T5|) . With small RR, some single par- 
ticle chains are also present, although shorter and less 
numerous. The sensitivity of plasticity index A to the 
rolling friction is likely to be explained by different rup- 
ture mechanisms, the importance of folding rearrange- 
ments growing with the level of rolling resistance. 



VIII. CONCLUSION 

To summarize, we have used numerical simula- 
tions to observe and characterize, at the macroscopic 
and microstructural levels, the consolidation behavior, 
in isotropic compression, of model cohesive powders. 
Macroscopic constitutives laws for quasistatic loading, 
unloading and elastic responses were shown to be reason- 
ably well approached. The material behavior was inves- 
tigated for a range of densities that is wider than in most 
simulation studies of cohesive granular materials. The 
consolidation process goes through three stages. In a first 
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regime, which is sensitive to the assembhng procedure, no 
plastic collapse occurs, as the agitation in the assembling 
process has stabilized a strong enough microstructure to 
withstand a finite pressure increase. The normal force 
distribution widens until a significant fi'action of contacts 
are on the verge of tensile rupture. The initial system ge- 
ometry, which changes very little in regime I, is that of a 
dense assembly of fractal blobs, with dimension dp taking 
the universal value associated with the aggregation pro- 
cess (here, ballistic) implemented in the sample prepara- 
tion stage. The blob size ^ (at most, between 5 and 10 
grain diameters in the present case) can be identified on 
studying density correlations. The subsequent consolida- 
tion behavior is remarkably independent on initial con- 
ditions, which merely determine where the intrinsic con- 
solidation curve in the <i>-P* plane is first met. The same 
curve is then followed whatever the initial conditions as 
the material is further compressed. This behavior cor- 
responds, at the microscopic level, to a gradual change 
of the blob size. The curve in regime II has the same 
shape as reported in the soil mechanics literature, and 
the consolidation pressure is a plastic threshold below 
which the material response is approximately elastic (like 
the behavior of a cohesionless granular material under 
isotropic load). Elastic moduli increase rapidly with con- 
solidation pressure or density. The cases of small RR or 
without RR should be distinguished from the situation of 
strong rolling resistance, although, in both cases, the mi- 
crostructure of loose packings might be viewed as denser, 
better connected regions joined by thin arms. In the first 
case, loose packings collapse when the tensile strength 
of contacts is overcome by the externally imposed forces, 
preferentially within the fragile junctions between adja- 
cent denser blobs. Systems with strong RR, on the other 
hand, contain single particle strands, which tend to fold 
without breaking in plastic compaction. While small RR 
systems gain very few contacts in the consolidation pro- 
cess, the coordination number might increase from nearly 
2 to 3 with large RR. Eventually, the material approaches 
a limiting, maximum density (regime III), as the pack- 
ing structure resembles that of a cohesionless system, for 
P* ^ 1 (albeit, typically, somewhat looser). The ab- 
sence of a similar upper limit of the density of cohesive 
packings in experiments for large P* is due to plastic 
deformation of contacts. 

The fractal blob size £,, depending on solid fraction $, 
is a central microstructural feature, based on which some 
scaling laws for elastic properties can be attempted. It 
is also tempting, beyond the qualitative description of 
the microstructural changes associated with the consoli- 
dation process, to try to predict the consolidation curve 
from such geometric data. Yet scaling laws only apply to 
a restricted part of the consolidation pressure interval. 

Our results, in many respects, emphasize important 
qualitative differences between cohesive and cohesionless 
granular assemblies. The existence of stable loose struc- 
tures and the consolidation phenomenon are the most 
important differences in macroscopic behavior brought 



about by cohesion. At the microstructural level, unlike 
in cohesionless packings, the typical values of intergran- 
ular forces, or the force distribution, are not as simply 
estimated in cohesive systems, in which attractive and 
repulsive contact forces of the order of tensile strength 
-Fo tend to compensate under low pressure. In particu- 
lar, compression cycles stabilize self-balanced force net- 
works with large compression forces. Unlike in granu- 
lar packings devoid of cohesion, the coordination number 
does not appear to be a significant state variable in co- 
hesive systems with low RR, as it hardly changes along 
the consolidation curve. With large rolling resistance, it 
witnesses, however, the formation of loops under com- 
pression. While cohesionless assemblies with low coordi- 
nation number usually contain many rattlers, all parti- 
cles in cohesive packings are connected to the same con- 
tact structure, which is rigid, but comprises lots of "dead 
ends" or "side arms", which might bear self-balanced 
forces but do not participate in the transmission of exter- 
nal stresses. Some of these new features can be summed 
up on remarking that loose powders are similar to gels 
as much as to granular packings with no cohesion. 



Our investigations should be pursued in several direc- 
tions. On the theoretical side, the connections between 
macroscopic properties and microstructure could be stud- 
ied more quantitatively. The behavior of loose cohesive 
packings under general stress states should be investi- 
gated. Thus one may determine whether such constitu- 
tive laws as the Cam-clay model [1| apply to the simu- 
lated material. And finally, more quantitative agreement 
with experiments and real materials should be sought. 
In spite of some obvious steps (e.g., one should simulate 
3D systems), this latter objective looks daunting. One 
major difficulty is the importance of hydrodynamic ef- 
fects at the assembling stage, when the microstructure 
and the fractal dimension of aggregates are determined. 
While we have bypassed this problem on implementing 
ballistic aggregation, it is necessary to investigate the 
behavior other possible kinds of aggregates, by dealing 
with some tractable model for hydrodynamic forces. It 
is hopefully possible to introduce some mechanics and in- 
tergranular interactions within the models used with ge- 
ometric aggregation rules (such as, e.g., diffusion-limited 
cluster-cluster aggregation). Then, another difficulty is 
that many parameters associated with the contact law 
(such as friction coefficient, rolling friction, rolling stiff- 
ness constant) should be identified for a real material to 
be investigated at the grain level. In this respect, the 
recent progress of experimental methods of microscopic 
investigation seems quite promising, as formerly inacces- 
sible parameters ruling interparticle contact mechanics 
are now beginning to be measured in model materials, 
thanks to particle-scale observation and micromanipula- 
tion techniques [13, HH, [H, [lE lli ■ 
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